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This  research  project  covered  the  period  1 October  1973  - 30  September 
1975.  The  fundamental  problem  addressed  was  the  discr  mination  between  earth- 
quakes  and  explosions.  At  the  outset  it  became  clear  that  data  of  high  quality 
would  be  needed  for  a quantitative  method  of  discrimination  cased  on  the 
seismic  source  mechanism.  Consequently,  a comparison  was  made  among  available, 
digitally  ’ecorded  seismographs.  Very  preliminary  results  are  presented  in 
Appendix  A It  appeals  that  the  _RG  instruments  are  definitely  superior  to  tne 
HGLP  instruments  although  both  sutler  from  unexplained  non-linearities.  Finite 
loop  gain  may  be  the  explanation  tor  the  GKO  instrument.  At  long  periods  th> 
fed-back  LaCoste-Romf erg  gravimeter  is  the  Lest  instrument  tested. 

To  retrieve  the  seismic  ..ource  mechanism  from  observed  spectra  one 
uses  the  concepts  ot  matched  filtering  and  deconvolution.  Before  the  beginning 
of  this  project  it  was  known  tnat  the  relationship  between  the  seismic  moment 
tensor  and  observed  spectra  is  a linear  one  (Gilbert,  1J71).  It  was  also 
known,  in  principle,  how  to  retrieve  the  moment  tensoi  (Gilbert,  1973).  The 
basic  ideas  were  refined  and  used  to  retrieve  the  moment  tensors  of  two,  larg-- 
deep  earthquakes  (Dziewonski  and  Gilbert,  1974;  Gilbert  and  Dziewonski,  1975). 

The  major  drawback  to  the  method  was  its  requirement  of  a large,  dense,  global 
network  of  stations.  Clearly,  a method  for  a regional  or  local  array  was  needed. 
Such  a method  was  discovered  as  the  project  drew  to  a close  and  t he  theoretical 
basis  for  it  is  presented  in  Appendix  B . In  theory,  one  needs  only  a single, 
horizontally  polarized  instrument  or  two,  vertically  polarized  instruments. 
Numerical  experiments  with  synthetic  data  indicate  that  as  few  as  5-10  verticals 


A practical  application  of  this  novel  method  is  presented  in  Appendix 
C where  it  is  shown  how  one  can  use  only  10  WWSSN  stations  to  retrieve  the 
mechanism  of  a deep  earthquake.  Digital  data  of  sufficient  quality  -tie 
instruments  must  be  well  calibrated  - were  simply  not  available  for  explosions. 
Therefore,  a practical  evaluation  of  this  method  awaits  the  forthcoming  SRQ 
data  which  is  expected  to  be  of  very  high  quality  (however,  see  Appendix  A). 

Fundamental  to  the  application  of  the  matched  filtering  method  is  the 
facility  to  calculate  very  accurate  synthetic  seismograms,  including  the 
effects  of  dissipation,  for  a broad  range  of  frequencies  and  wave  numbers. 

Among  all  known  methods  the  classical  procedure  of  summation  of  normal  modes 
was  found  to  be  the  most  reliable.  To  compute  the  required  normal  mode 
eigenfrequencics  and  eigenfunctions  is  not  a trivial  task,  even  after  20  years 
of  theoretical  and  numerical  effort.  We  have  found  that  the  classical 
Rayleigh-Ritz  procedure  is  the  cheapest  and  most  accurate.  All  eigendata  with 
periods  between  40  sec  and  1 hour  - about  5000  modes  - can  be  computed  for  a 
few  thousand  dollars.  The  present  programs  are  designed  to  produce  the  com- 
plete spectrum  for  all  periods  greater  than  20  sec.  The  theoretical  basis  for 

this  variational  calculation  is  given  in  Appendix  D. 

At  periods  on  the  order  of  a few  tens  of  seconds  the  concept  of  standing 

waves,  valuable  at  longer  periods,  is  better  replaced  by  the  concept  of 

traveling  waves.  An  exact  theoretical  traveling  wave  representation  is  pre- 

4 

sented  in  Appendix  L.  Such  a representation  is  very  desirable  for  applying 
the  matched  filtering  methods  to  a regional  or  local  network. 

In  summary,  a novel  method,  l>as:d  on  matched  filtering,  has  been  found 
for  retrieving  the  seismic  source  mechanism,  the  moment  tensor,  from  a sparse 

The  moment  tensor  can  be  examined  for  its  deviatoric 


network  of  instruments. 
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(e^thquake-like ) and  isotropic  (explosion-like ) components.  Almost  all 
earthquakes  are  thought  to  have  a nearly  completely  deviatoric  moment  tensor, 
and  almost  all  explosions  to  have  an  isotropic  one.  Thus,  the  unique  partition 
seismic  moment  tensor  into  its  deviatoric  and  isotropic  parts  provides 
a quantitative,  unambiguous  method  for  discriminating  between  earthquakes  and 
explosions.  The  matched  filtering  method  for  retrieving  the  seismic  moment 
tensor  works.  It  is  new  and  its  limitations  have  not  been  explored.  Clearly, 
it  is  scientifically  most  desirable  that  the  full  potential  of  this  method 


be  assessed. 


Appendix  A 

We  present  some  preliminary  comparisons  among  instruments  for 
the  WWSSN , HGLP , LRO  and  IDA  networks.  The  comparisons  are  incomplete 
yet  adequate  enough  to  show  some  of  the  features  of  the  different  systems. 

In  Figure  A1  we  have  four  seismograms  of  the  Korean  earthquake 
(m^  = 6.5)  on  September  29,  1973.  The  four  instruments  (all  verticals) 
are:  the  Goodkind-Prothero  superconducting  gravimeter,  the  La  Coste- 

Romberg  gravimeter  (the  IDA  instrument),  the  HGLP  vertical  and  the  WWSSN 
long-period  vertical  (recorded  digitally).  The  HGLF  vertical  is  located  in 
Albuquerque,  N.M.  and  the  other  three  at  the  Pincn  Plat  Observatory,  Calit. 

The  two  gravimeters  are  band-pass  filtered  to  enhance  (xlCO)  the  acceleration  in 
the  band  1 min-30  min.  The  HGLP  vertical  has  the  standard  "noise  notch" 
filter  and  the  WWSSN  vertical  has  the  standard  response.  Rayleigh  wave 
packets  through  R6  are  clearly  visible  on  the  two  gravimeters. 

In  Figure  A2  we  have  the  spectra  for  the  first  day  after  the 
earthquake  from  the  seismograms  in  Figure  A1  (the  WWSSN  spectrum  is 
from  a hand-digitized  record  of  the  Albuquerque  WWSSN  instrument).  The 
filters  on  the  two  gravimeters  are  virtually  identical  in  the  pass  band 
and  so  are  the  spectra  of  the  two  instruments.  The  similarity  between 
the  two  of  the  spectral  peaks  and  troughs  indicates  a S/N  ratio  approaching 
50  db.  At  low  frequencies , near  10  cph,  the  WWSSN  vertical  has  a S/N 
ratio  of  at  least  25  db,  equal  to,  if  not  better  than,  the  HGLP  instrument. 

This  comparison  shows  that  a digitally  recorded  WWSSN  vertical  is  at 
least  marginally  preferable  to  the  "A"  channel  of  the  HGLP  vert-cal  and 
that  the  two  gravimeters  are  some  25  db  better  than  either  of  the 
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other  two  instruments.  It  is  not  clear  whether  the  other  two  instruments 
respond  to  ground  noise  at  lor.g  periods  because  the  dynamic  range  of  the 
recording  system  is  taken  by  the  large  amplitude,  short  period  motion. 

In  the  present  configuration  ground  noise  could  be  as  much  as  25  db  below 
instrument  noise  at  long  periods  for  both  i he  WWSSN  and  HGLF  instruments. 

In  Figure  A3  we  have  spectra  for  the  second  day  of  the  earthquake. 
The  WWSSN  spectrum  has  been  omitted.  A glance  at  Figure  A1  shows  that 
the  ambient  noise  level  is  reached  sore  6 hours  after  the  earthquake 
for  this  instrument . A number  of  long  period  fundamental  modes  and  high-Q 
overtones  are  clearly  present  in  the  spectra  of  the  two  gravimeters  and, 
to  a lesser  extent  , in  the  HGLP  spectrum.  At  10  cph  the  spectral  level 
has  dropped  15  db  from  day  1 to  day  2 for  the  gravimeters  and  30  db  for 
the  HGLP  instrument.  This  demonstrates  the  nonlinearity  of  the  HGLP 
instrument.  Large  amplitude,  short  period  motions  are  non-linearly 
"aliased"  to  long  periods.  As  ."he  short  period  motions  decay,  and  they 
dec^y  more  rapidly  than  the  long  period  motions,  the  aliased  long  period 
motions  also  decay  more  rapidly  than  the  true  long  period  motions.  Thus, 
much  of  the  long  period  spectral  energy  in  the  HGLP  spectra  in  Figures 
A2  and  A3  is  non-linearly  aliased  short  period  spectral  energy.  This 
intolerable  situation  makes  the  HGLP  instrument  unacceptable  for  long 
period  studies. 

To  improve  the  comparison  of  the  HGLP  and  La  Coste  instruments 
at  long  periods  we  turn  our  attention  to  the  boom  channel,  channel  "B", 
of  the  HGLP  instrument. 

In  Figure  A4  we  have  the  tide  channel  (flat  in  acceleration)  of 
the  La  Coste  gravimeter  operating  in  Nana,  Peru  and  in  Figure  A5  the 
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"B"  channel  of  the  HGLP  instrument  at  l.a  Paz,  Bolivia,  both  for  the 
Solomon  Islands  earthquake,  h = 7.8,  July  20,  1975.  Notice  that  the 
La  Coste  is  recorded  as  an  acceleration  indicator  and  the  HGLP  as  a 
displacement  indicator , so  that  the  two  records  have  tides  that  are  cf 
opposite  signs.  In  both  instruments  the  large  tidal  signal  takes  most 
of  the  dynamic  range  of  the  recording  system,  leaving  a few  tens  of 
least  count'  for  this  rather  large  event.  The  desirability  of  band- 
pass filtering  is  obvious,  and  in  figure  AG  we  show  the  La  Coste  record 
amplified  <100  in  the  band  1 min  - JO  min.  This  ampl i i icat ion  is  done 
by  active  filters  bej eve  digitization  so  '.hat  the  earthquake  has  a few 
thousands  of  "least  count",  or  about  70  db,  dynamic  range.  Spectra 
corresponding  to  Tigures  A5  and  A6  are  shown  in  figures  A7  and  A8. 

Two  spectra  are  presented  in  each  figure, one  for  the  24  hr  period  before 
the  earthquake  and  the  other  for  the  first  24  hours  after  the  earthquare. 
At  long  periods  the  La  Coste  gravimeter  has  a S/N  ratio  for  this  event 
of  about  50  db  while  the  HGLP  S/N  ratio  is  nowhere  more  than  15  di . The 
noise  level  must  thus  be  at  least  35  db  abo*  e.  ground  noise;  some  of  this 
is  presumably  least-count  noise,  but  nonline irity  in  the  sensor  (folding 
higher-frequency  ground  noise  to  lower  frequencies)  and  noise  in  the 
instrument  are  also  contributing. 

Data  for  the  SRO  instruments  has  only  very  recently  become 
available.  Therefore,  a detailed  comparison  of  the  SRO  instrument  with 
others  is  not  now  possible.  We  have  in  figure  A9  the  seismogram  of  the 
Solomon  Islands  event  recorded  on  the  SRO  vertical  at  Albuquerque. 
Although  this  instrument  is  fedback  (electromagnetically ) the  output 


is  filtered  with  a sharp  rolloff  at  long  periods  before  digitization. 

We  have  low  passed  the  original  record  with  a corner  at  100  sec.  to 
obtain  Figure  A9.  Spectra  for  the  day  before  and  for  the  first  24  hours 
of  the  earthquake  are  shown  in  Figure  MO.  The  S/N  ratio  at  long 
periods  is  40  to  50  db,  about  equivalent  to  the  La  Coste  gravimeter. 

The  lack  of  structure  in  the  earthquake  spectrum  for  periods 
longer  than  about  300  sec  is  attributable  to  the  unusually  large  Rayleigh 
wave,  Rl,  for  the  main  shock.  This  signal  represents  an  instrumental 
problem;  whether  the  active  filters  are  saturated  or  there  are  difficulties 
in  the  gain  ranging,  Rl  as  recorded  on  the  SR0  instrument  appears  to 
be  anomalous.  Excising  Rl  from  the  seismogram  leads  to  the  spectra 
shown  in  Figure  All  (the  noise  spectrum  is  the  same  as  in  Figure  A10). 

The  spurious  long-period  noise  disappears,  and  moo»  peaks  can  be  seen 

out  to  470  seconds.  Saturation  is  a problem  with  any  seismographic 
instrument.  For  the  La  Coste  gravimeter  saturation  at  Rl  appears  to  be 
a problem  at  teleseismic  distances  for  magnitudes  M > 7lj  . For  the  SR0 
instrument  the  problem  of  saturation  is  not  yet  well  understood. 

At  the  longest  periods  at  which  modes  are  visible  on  the  SRO 
spectrum,  the  S/.N  ratio  is  only  25  db,  compared  with  40  db  for  the 
La  Coste.  Indeed,  the  S/N  ratio  for  the  SRO  instrument  begins  to  deteriorat 
for  periods  longer  than  300  seconds,  whereas  the  La  Coste  response  remains 
good  out  to  600  seconds.  As  a result,  the  longest-period  mode  unequivocally 
visible  on  the  SRC  spectrum  is  g^l3  at  474  seconds;  c the  La  Coste 
spectrum  ^S^  is  visible,  at  a period  of  1190  seconds.  Another  indication 

of  the  quality  of  the  La  Coste  gravimeter  is  presented  in  Figure  A12. 
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Wa  have  four  sections  of  tha  seismogram  (amplified  *100  in  the  band 
1 min  - 30  min)  for  the  Solomon  Islands  event  spaced  roughly  half  a 
day  apart.  The  semi-diurnal  tide  is  evident  in  this  figure.  The 
start  times  of  the  records  are  -8.5,  tl6,  *28,  and  +41  hours  with 
respect  to  the  origin  time  of  the  earthquake.  Up  to  2 days  after  this 
M-7.8  event  there  is  still  evidence  of  the  free  oscillations  of  the 
Eavth. 

In  conclusion,  it  appears  that  the  La  Coste-Ronberg  gravimeter 
is  the  instrument  of  choice  for  a sparse,  long-period  network.  If 
Project  IDA,  or  a similar  project  also  using  La  Coste  instruments,  is 
not  supported  then  seismologists  will  not  have  high-quality,  long-period, 
digital  data.  Obviously,  such  data  are  very  desirable.  The  only  ->ther 
instrument  that  might  possibly  be  useful  for  long-period  seismology  is 
the  SRO  instrument  (Peterson,  et  at  . , USGS  preprints,  July-August  1375). 
At  present  little  is  known  about  the  noise  characteristics 
of  the  SRO  instrument  at  very  long  periods,  nor  are  the  nonlinearities 
caused  by  large  signals  well  understood.  It  thus  seems  that  the  fedback 
La  Coste  is  the  only  proven  low-noise  vertical  seismometer  at  very 


long  periods. 
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FIGURE  CAPTIONS 
Append  iv.  A 

A1 . Seismograms  from  the  North  Korean  earthquake  of  September  29,  197 
recorded  on  a La  Coste  gravimeter  and  superconducting  gravimeter, 
and  standard  Press-Ewing,  at  Pif.on  Elat  Observatory,  and  ar  HGLP 
at  Albuquerque. 

A2.  Spectra  from  the  first  24  hours  after  the  North  Korean 
earthquake . 

A3.  Spectra  from  the  second  21*  hours  after  the  North  Korean 
earthquake . 

A4.  Solomon  Islands  earthquakes  of  July  20  and  21,  1375  recordec 
on  the  tide  channel  cf  the  La  Coste  gravimeter  at  Nana,  Peru. 

AS.  The  Solcnon  Islands  earthquakes  on  the  B Channel  of  the 
HGLP  at  La  Pan,  Bolivia. 

A6 . The  Solomon  Islands  earthquakes  on  the  Nani  La  Coste  gravimeter 
with  the  band  .3  to  10  mHz  amplified  by  100. 

A7.  Power  spectrum  of  the  Solomon  Islands  earthquakes  on  the  HGLP 
B channel.  Light  line  is  for  24  hours  before  the  earthquake;  heavy 
line  is  for  24  hours  after. 

A8.  Power  spectrum  from  the  NNA  x 100  record. 

A9.  Solomon  Islands  earthquakes,  recorded  on  the  SPG  vertical  at 
Albuquerque;  data  have  been  digitally  lopassed  to  remove  frequencies 


above  10  mHz. 
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A10.  Power  spectra  of  the  seismogram  in  A9.  (Light  line  is  spectrum 
of  previous  24  hours.) 

All.  power  spectrum  of  signal  in  A9,  but  with  R1  excluded, 

A12.  Signal  from  xlOO  output  of  NNA  La  Coste  at  times  before  and 

after  Solomon  Islands  earthquakes,  to  show  rate  of  signal  decay. 
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Summary 

In  theory,  a single  horizontally  polarized  seismometer  can  t< 
used  to  find  the  six  independent  elements  of  the  seismic  moment  tens, 
of  a buried  point  source,  provided  that  the  instrument  is  neither 
longitudenally  nor  transversely  polarized.  Also,  two  vertically 
polarized  seismometers  can  be  used,  provided  that  the  epicenter  does 
not  lie  on  the  great  circle  through  the  two  instruments.  These  resu; 
form  the  theoretical  basis  for  a procedure  for  retrieving  the  source 
median  ism  from  a sparse  seisrr.ographic  network. 


seismic 


Let  the  six  independent  elements  of  the  seismic  moment  rate 
spectrum  be  f (w ) = (t^(w)  , . . . , f^w))7  and  suppose  that  P 

spectra  (records)  u(u>)  = (u^w) UpC-))7  have  been  observed. 

The  relationship  between  u (w)  and  f_(^)  is  a linear  on«  (Gilbert,  1971, 
1973:  Dziewonski  and  Gilbert,  1974;  Gilbert  and  Dziewor.ski,  1975,  hereafter 
referenced  as  d ) 

u(w)  - H(w)  • f(u>) 

The  1x6  matrix  is  a functional  of  the  mechanical  structure 

of  the  Earth  and  can  be  regarded  as  the  spectral  ti  ester  mati ix  or 
system  function  that  relates  output  u(w)  to  input  f(w)  . Let  the 
p^-  row  of  H(ui)  be  h7p(w)  . 

The  six-vector  h (w)  can  be  written  as  the  sum  of  normal 

r 

modes  (M;  2.1.24,  2.1.28) 

Vw)  = £ Ck(w)  ^2) 

where  A specifies  the  excitation  and  amplitude  of  the  k— mode  for 
the  p— record  , is  the  resonance  function  of  the  k— mode,  and 

represents  the  effect  of  truncation  and  the  response  of  the  p— 
instrument.  Each  element  of  hp(w)  is  the  spectrum  of  a seismogram 
caused  by  a unit  element  o:  the  moment  rate  tensor,  a delta  function 
in  time.  The  observed  seismogram  at  the  p—  instrument  is  a linear 


& 
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combination  of  che  six  seismograms  h (u)  and  the  six  coefficients  in 

“P 

the  linear  combination  are  the  six  elemr  ,.s  of  f(u)  . 

Suppose  that  our  model  of  the  Earth  is  good  enough  to  permit 

us  to  ignore  the  difference  between  real  and  calculated  h (w)  . Then 

~P 

we  can  seek  to  solve  (1)  for  f_(w)  . At  low  frequencies  the  spectral 
peaks  in  b^p(w))  are  sufficiently  well  separated  to  cause  spectral  gaps, 
frequencies  where  there  is  little  or  no  information  abcut  f(w)  . 
However,  it  is  generally  believed  that  f_(ui)  is  a smooth  function  of 
uj  at  low  frequencies,  so  smooth  that  it  can  be  taken  constant  over 
a frequency  band  embracing  many  modes.  Therefore,  define  set  J of 
discrete  frequencies  um 

(3) 

U.J  - 6u.  s u>i  < Wj  + 6ui  i i = ij  , ij  + 1 i + I -1 

and  replace  f_(u>)  by  f_Coj j ) . There  are  now  I • P equations  for 
the  six-vector  f 


Ej=  «j  * <4> 

and  we  solve  (4)  by  applying  the  classical  method  of  least  squares 

s!j  • ’Sj  ■ b“j>  i-Tj  = h"  . H J (5) 


where  the  superscript  H denotes  hermitean  transpose. 


In  forming  (5)  we  cross -correlate  u with  each  element  of  h 

P -P 

H 

(multiply  by  h (<u)  ) . This  operation  is  conventionally  termed  matched 


filtering  and  is  an  operation  to  enhance  the  signal  being  sought.  The 





£ 3 


result  is  summed  over  the  frequencies  in  set  J to  give  (5).  In  order 


to  solve  (5)  for  f we  require  that  .**  have  rank-6  . Thus  I • P>  6 


is  a necessary  condition.  For  a dense  network,  P >>  1 , I , the  number 
of  discrete  frequencies  in  set  J , can  be  small.  Alternatively  if 
12  6 it  appears  that  we  can  have  P = 1 and  still  maintain  rank-6 
for  • To  explore  this  possibility  we  examine  the  eigenvalues  of 

• Without  loss  of  generality  we  take  R (dJ  = l 

kp 

Consider  a single,  vertically  polarized  accelerometer . In 
epicentral  spherical  coordinates  the  location  of  the  receiver  is 
(r  » 6 , ♦)  • An  inspection  of  (M,  2.1.30)  shows  that  the  six  vector 
(2)  for  vertical  polarization  (r-component ) can  be  written 


A = (*  • S)  u(r) 


(6) 


where  jt  is  a 6 x 4 matrix  whose  ncr.-zero  elements  are 


*11  = *22  = *32  = 1 ’ *23  = "*33  = cos2$  » *44  = cos* 


'44 


(7) 


♦ 54  = sin#  , = 2 sin2# 


63 


and  _S  is  a 4-vector  with  components  (M;  2.1.30) 


.0  v0 


«I  ■ s2  ■ 4 X°t  • S3  ■ *4  XJ  , 


2 .2 


2f^ 

4 l 


C8) 


I 


Substituting  (6)  into  (2)  gives, 

h(r  ,<*>)=£•  P(r  , w)  ; P(r  , u,)  = £ S,  C,  U)  U,  (r)  (9) 

k -k  k k 

where  P(r  , u ) is  a 4-vector.  For  ,yf  we  have 

zz  J 

£ T (10) 

where 

£ £*<r  1 (**!>  pV  . w.)  (11) 

“ i 

and  the  asterisk  denotes  complex  conjugate.  In  (10)  the  4x4  he^mitean 
matrix  & is  limited  to  rank-4  and,  therefore,  so  is  . Consequently, 

the  6x6  hermitean  matrix  is  singular,  and,  not  surprisingly,  the 

moment  rate  tensor  cannot  be  retrieved  from  the  spectrum  of  a single, 
vertically  polarized  accelerometer . However,  if  I , the  number  of 
discrete  frequencies  in  set  J , is  large  enough  (la  4 is  necessary) 
then  & can  have  rank-4.  We  shall  assume  that  Ji®  is  rank -4. 

Consider  two,  vertically  polarized  accelerometers  with  cooruinates 
(6j  > ana  (e^  , $2)  . In  an  obvious  notation  (10)  becomes 

= £l  ' — 1 *—  1+  — 2 (12) 

If  then  ^ and  (11)  becomes 

Jifj  = * • (£1  t # 2)  • * T (13) 

making  singular.  Also,  if  |<j>^  - $2|  = tt  , £ ^ and  £ are  the 

same  except  for  the  sign  of  column  4.  If  we  change  the  sign  of  column  4 
and  row  4 of  djPj  we  have  (13)  again.  Therefore,  if  the  epicenter  lies 


on  the  great  circle  through  the  two  vertical  Instruments,  is 

singular.  Included  here  is  the  special  case  0^  = 0 , tt  or  = 0 , it 
By  a proper  choice  of  coordinates  we  can  always  have 
~ *2  = *1  ~ * in  -12)*  Assuming  that  & 1 and  a?2  both  have  rank-4 
we  deal  with  the  matrix  *_  (*)  • * T(*)  + *_  (-*)  . $ T(.$)  which 

has  rank-6  unless  * = 0 , tt/2  , it  . Therefore,  if  the  epicenter 

does  not  lie  on  the  great  circle  through  the  two  instruments,  .Yf  is 

* IT  J 

non-singular  and  (5)  can  be  solved  for  f(oj  ). 

J 

This  result  is  important  because  it  shows  that  a sparse  global 
network  of  vertically  polarized  instruments  can  be  used  to  retrieve 
the  seismic  moment  rate  tensor.  Buland  and  Gilbert  (1976)  have  shown 
that  using  ten  WWSSN  stations  leads  to  a satisfactory  result  for  rr,^  ~ 7 
To  consider  horizontally  polarized  instruments  we  must  take  into 
account  toroidal  as  well  as  spheroidal  modes.  We  introduce  the  6x2 
matrix  whose  non-zero  elements  are 

*21  = "f31  = sin2*  * *61  = -2cos2$  , *42  = sin*  , = -cos* 

and  the  2-vector  T (M;  2.1.31) 


X 


l 

l 


In  terms  of  f and  T , A 


(14) 


(15) 


in  (M;  2.1.28)  for  toroidal  modes  is 


A *-8  csce  3^  (f  • T)W(r)  ♦ *_  30  (f  • T)  W(r)  (16) 

and 

h^(r  , ui)  = -3  , • 53  csc  6 T^C^(u) )W^(r ) = • Q.  esc  0 

_ 9 — k ~ 

(17) 

h^t  ,<*>)=»•  53  30  TkCk(u1)Wk(r)  = I • Q ' 

~ k — 

hj  and  for  spheroidal  modes  are 

^(r.w)  = • * 3e  D = i • D'  ; D(r  , w)  = 53  ^ Ck(w>  V (r) 

k 

118) 

hg(r_  , u)  = 3^  ♦ • D_  esc  0 = ♦_"  • D esc  0 
Thus  the  complete  "synthetic  seismograms"  are 

h_2(r_,  oj)  = «2  • R2(u)  , £2  = * « -f_'  » R 2 = D*  ® £ csc  0 

(19) 

hg(r_ , u)  = C w ) . = ♦'  « i 1 Rg  = CSC  0 D • q/ 

where  the  6x6  matrices  fi.2  and  are  functions  only  of  $ and  the 
6-vectors  R_2  and  are  functions  of  r,  0 and  u.  For  both  fi2  and 
the  fifth  column  is  proportional  to  the  third  and  the  sixth  to  the  fourth. 
Also,  columns  1 and  2 of  n,  are  zero.  Thus  n„  has  rank -4  and 

Z?  —2 

flg  has  rank-2.  When  we  sum  over  set  J to  obtain 
Jfj  we  use  the  6x6  matrices 


ft7 


z 


R^(u.)  R^u.)  ; g , Y = 2 


, 3 


(20) 


H 

and  we  have  We  assume  ,jg  to  have  rai.k-6.  Although 

this  assumption  will  be  supported  for  w-bands  that  include  multiplets 
for  several  values  of  £ approaches  singularity  for  large  £ . 


This  is  a result  of  ||  D*I|/||D_I!  = 0 (£)  for  large  £ . The  sane  is 
true  for  Q_  . This  means  that  ^q2  aPProac^es  rank-4  as  an  upper  left 
4x4  block,  approaches  rank-2  as  a lower  right  2x  2 block, and 

approaches  rank-2  as  an  upper  right  4x2  block.  Physically,  this 
decomposition  is  a result  of  spheroidal  modes  dominating  the  ©-component 
and  of  toroidal  modes  dominating  the  ^-component  for  large  £ . 

We  now  consider  a single,  horizontally  polarized  instrument 


oriented  at  an  angle  a with  respect  to  the  6 -vector.  In  terms  of 
(19)  and  (20)  .X?Tis 


.*j  = cos2o  n2  • 3^2  • n7  + cosa  sina(n_2  • .^23  • nj  + 


+ -%2  * & + sin2a  Sg  *^3  * Sj 


(21) 


In  general,. will  be  non-singular.  However,  if  a = 
polarization  , or  -it/ 2 , transverse  polarization, 

— J 

because  n2  and  ^ are  singular.  Also,  the  matrices 


: 0 , longitudenal 
will  be  singular 


ey 


have 


rank-2  for  6 = 0 , n . This  means  that  a source  directly  beneath 
the  receiver  or  its  antipode  cannot  be  retrieved.  Otherwise,  the 
moment  rate  tensor  can  be  retrieved  from  a single  , horizontally 
polarized  instrument.  As  in  the  previous  example,  for  two  vertical 
instruments,  it  is  necessary  to  sum  over  an  w-band  containing  multiplets 
for  several  values  of  t , in  order  that  •#  have  full  rank,  and  it  is 

— By 

assumed  that  f(w)  is  nearly  constant  in  each  w-band.  This  result 
remains  true  for  large  i even  though  becomes  singular.  Since 

&22  becomes  an  upper  left  4x4  block  we  can  replace  by  £ in 
(21).  Similarly,  we  can  replace  ^ by  L • Let 

2 2 

= cos  + sin  <*.#33  + cosa  sina 

(22) 


U = * • V 


For  large  l (21)  becomes 


Jffr  2.  • 4*  • nT  (24) 

=*I  — — _ 

In  (24)  we  assume  that  I , the  number  of  discrete  frequencies  in  set  J , 
is  large  enough  to  make  i#  have  rank-6  . Since  det  £)_  * 0 (actually 

det  R = 8)  we  see  that  Jtfj  has  rank-6  . Thus,  even  at  short  periods, 

the  moment  rate  tensor  can  be  retrieved  from  a single  horizontal 

instrument  unless  a=0  , 1/ 2 op  0=0  ,n  , 

In  practice,  a seismographic  station  has  two  horizontal  instruments, 

in  which  case  it  is  clear  that  the  moment  rate  tensor  can  be  retrieved. 


Moreover,  a standard  installation,  consisting  of  one  vertical  and  two 
horizontal  instruments  certainly  enables  the  retrieval  of  the  moment 
rate  tensor.  Here,  the  only  exclusion  is  3 = 0 , it  . 

The  foregoing  examples  demonstrate  theoretically  that,  except 
in  special  circumstances,  the  moment  rate  tensor  of  a buried  point 
source  can  be  retrieved  from  the  spectra  of  two  vertical  accelerometers 
or  from  the  spectrum  of  one  horizontal  accelerometer.  From  these 

theoretical  results  we  can  easily  infer  that  a network  of  a small 
number  of  instruments  can  be  used  to  retrieve  source  mechanisms  on 
a routine  basis.  The  ability  to  achieve  such  retrievals  makes  possible 
some  interesting  research  projects. 

The  method  presented  here  is  an  extension  of  the  concept  of 
matched  filtering  (see,  for  example,  Robinson,  1967,  pp . 259-264). 

The  matched  filters,  h(r  , u >)  in  (2),  are  the  best  linear  filters 
in  that  they  maximize  the  signal/noise  ratio.  For  Gaussian  noise 
they  are  optimum. 

Although  we  have  obtained  h_(r  , u)  in  (2)  by  summing 
normal  mode  multip.lets,  it  should  be  emphasized  that  the  method  of 
retrieval  is  independent  of  the  procedure  used  to  obtain  h_(r  , w)  . 

Any  procedure  for  generating  synthetic  seismograms  can  be  used  to 
obtain  h_(r  , <d)  . Therefore,  matched  filtering  for  the  seismic 
moment  tensor  can  be  done  globally,  regionally  or  locally,  depending  on 
the  magnitude  of  the  seismic  source  and  the  configuration  of  the 


network. 
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ABSTRACT 

By  a process  of  matched  filtering  it  is  possible  to  deconvolve 
a small  number  of  acceleration  or  strain  records  for  the  moment  rate 
tensor  of  a seismic  source.  Theoretically,  one  horizontal  accelerometer 
(or  strain  meter)  or  two  vertical  accelerometers  is  sufficient.  Practically, 
five  to  ten  records  can  be  shown  to  suffice.  Specifically,  the  mechanism 
of  the  Colombian  event  of  July  31,  197C  can  be  retrieved  from  a sparse 
network  of  ten  WWSSN  vertical  instruments.  With  currently  available 
instrumentation  it  should  be  possible  to  discriminate  unambiguously 
between  an  earthquake  and  an  explosion  larger  than  magnitude  6 at 
teleseismic  distances.  Improved,  digitally  recorded  networks 
lead  to  a decrease  in  the  threshold  of  discrimination. 


can 


We  consider  the  moment  rate  tensor  for  a buried  point  source. 
The  problem  of  retrieving  the  elements  of  the  moment  rate  tensor 
(Gilbert  and  Dziewonski,  1975,  Section  2.3;  hereafter  referenced  as  M) 
has  been  discussed  by  Dziewonski  and  Gilbert  (1974).  Their  method,  an 
extension  of  a suggestion  of  Gilbert  (1973),  depends  on  the  orthogonal! 
relations  for  the  spherical  harmonics,  and  implies  a dense  global  array 
of  receivers.  The  paucity  of  first  quality  long  period  instruments 
motivates  a search  for  a method  requiring  only  a sparse 
network. 

Following  (M;  2.1.17),  let  the  Fourier  transform  of  the 

moment  rate  tensor,  M^)  , be  written  as  the  six-vector,  f(fc)  . 

Then  the  Fourier  transform  of  the  p—  displacement  (or  strain,  etc.) 

record,  due  to  f(w)  , at  the  i—  frequency,  U (w.)  , mav  be  written 

P i 

as  a sum  of  normal  nodes  (M;  2.1.24,  2.1.28). 


VV  = ? 4jp  • l <“i>  ck<“i>  VV 


(1) 


Where  specifies  the  excitation  and  amplitude  of  the  k—  mode 


th  . 


instrument;  is  the  resonance  function  of  the  k^— 


at  the 


normal  mode;  and  is  the  instrument  response  and  truncation  effect 

at  the  p—  instrument.  Define  matched  filters 


h («.)  = T A,  C.  (u.  ) R (u.)  . 

~P  1 v K 1 P 1 


For  P records  we  have  P equations  in  six  unknowns;  f , . . . , f 


U(uk)  = H(u^)  • £(u.) 


where  the  p—  element  cf  U(u..)  is  U (w.)  and  the  p^-  row 
r — l p i ^ 

T 

of  H(u.)  is  h (u.)  • Clearly,  if  P > 6 we  can  solve  (3)  by 

ZZ  1 1 J 

least  squares  for  ftw^)  , the  moment  rate  tenser  at  frequency 
point  . 

Because  of  ground  and  instrumental  noise,  and  splitting  and 
uncertainties  in  Q structure,  it  is  unlikely  that  £(u^)  will  be  well 
determined  for  a small  number  cf  stations.  However,  from  M( Figures  6 and  7) 
we  suspect  that  f is  a smooth  function  of  frequency.  Therefore,  define 
set  J of  discrete  frequencies  such  that  Uj  - 6u  s u.  s Uj  + 6u  , 

i = l,2,...,I.  Assume  that  f_(uij)  = f(w^ ) = . . . = f_(u>j)  . 

We  now  have  P • I equations  in  six  unknowns 


U = H • f(u.)  . 
— d rd  — J 


As  I , the  number  of  frequency  points  in  set  J , may  be  large 
there  appears  to  be  no  reason  that  P cannot  be  unity.  To  explore 


C 3 


the  limits  of  the  method  we  have  examined  the  eigenvalues  of  the  normal 
matrix 

•Tj  = " * ' LLj  (s) 

where  the  superscript  * means  hermitean  transpose.  A detailed 
analysis  is  beyond  the  scope  of  this  letter  and  will  appear 
elsewhere.  The  results  will  be  stated  without  proof.  Two  vertical 
accelerometers  (not  lying  on  a great  circle  through  the  epicenter)  or 
one  horizontal  accelerometer  (not  long: tucinaily  or  transversely 
polarized)  give  full  rank  to  . Practically,  synthetic  numerical 

experiments  indicate  that  a minimum  of  five  to  ten  stations  is 
satisfactory. 

For  large  P we  examine  the  relationship  between  the  above  method 
and  the  method  described  in  (M;2.3).  for  simplicity  assume  P (w)  = 1 . 
Multiplying  equation  (4)  by  H (explicity  forming  the  normal  equations) 
and  expanding  the  sums 


2 {s  [ ? "‘'“dj  [E  4 w]}  • 


C6) 

,)  • 


Rearrange  equation  (6) 


& [S  <(“d  W]  b.P 

*{fj  f?  ['Z  bp  4]}  • 


C7) 


and  define 


ut+6uj 

w >. 


7 is  _ 

Ck(w)  V (u)du  = £ CK(u)i)  v (u)i)  • 


w t-6oj 
J 


Using  the  orthogonality  relation  (M ; 2.3.2)  among  modes  (valid  for 
a dense  array)  rewrite  equation  (7)  as 


E MU  > A, 

k P 


= I E 1 A.  aJ  1 . fi 

lp,K  ^ V'  ‘ 


f(0)j)  . 


Equation  (9)  differs  from  M;  2.3.2-2.3.10  only  in  the  definition  of 
I(v)  . However,  equation  (8)  shares  with  its  counterpart  (M;  2.3.7) 
the  insensitivity  to  Q gained  by  integrating  across  the  resonance 
function.  The  fundamental  difference  between  matched  filtering 
(equation  (6)  ) and  stacking  (equation  (9)  and  (M;2.3.8))  is  that  matched 
filtering  includes  all  the  cross  terms  among  modes  (does  not  depend 
on  the  orthogonality  relation). 

The  method  presented  here  employs  an  extension  of  the  concept 
of  matched  filtering  (see,  for  example,  Robinson,  1967,  pp.  259-264). 

The  matched  filters,  ) , in  (2),  are  the  best  linear  filters 

in  that  they  maximize  the  signal/noise  ratio.  Intuitively,  we  remove 
noise  and  scattered  energy  from  a record  by  forcing  our  prejudice  on 
each  record  that  it  be  a linear  combination  of  six  synthetic  records. 

The  data  are  then  frequency  averaged  over  set  J (sum  on  i in  equation 
(6)  ) to  reduce  sensitivity  to  splitting  and  Q structure  and  to  stabilize 
•%j  . Spatial  averaging  (sum  on  p in  equation  (6))  further  stabilizes 
the  system.  Finally,  the  P records  are  deconvolved  for  f(<uT)  by 
inverting  . 

“J 


C.5 


As  an  example  of  the  method  we  have  plotted  (figure  1)  the  raw 
M component  derived  by  matched  filtering  from  a ten  record  subset 
(vertical  first  day  data  from  stations  GUA,  JER,  KIP,  CHG,  GIE,  GDH, 

PTO,  NAI,  NAT,  KBL)  of  the  Colombian  data  set  (M;  Table  1)  and  M 
redrawn  from  (M,  Figure  6).  The  agreement  is  apparent  and,  considering 
the  quality  of  WWSSN  data  at  very  long  periods,  acceptable.  Figure  1 
emphasizes  the  demand  for  quality  data  if  P is  to  be  small. 

Given  the  availability  of  quality,  long  period,  digital  data, 
matched  filtering  provides  a method  for  routinely  and  rapidly  determining 
the  complete  long  period  source  function  of  any  event  larger  than  about 
magnitude  6 at  teleseismic  distances.  This  should  greatly  facilitate 
studies  of  stress-release  mechanisms.  Furthermore,  the  implications 
for  seismic  discrimination  are  promising.  Having  calculated  the  complete 
long  period  source  function  of  an  event,  the  discrimination  between 
an  earthquake  and  an  explosion  is  entirely  unambiguous. 


Although  we  have  obtained  h^( ) by  summing  normal  mode 
multiplets,  it  should  be  emph  :ized  that  the  method  of  retrieval 
is  independent  of  the  procedure  used  to  obtain  the  matched  filters. 
Any  procedure  for  generating  synthetic  seismograms  can  be  used. 


Therefore,  matched  filtering  for  the  seismic  moment  tensor  can  be 


done  globally,  regionally  or  locally,  depending  on  the  magnitude  of 
the  seismic  source  and  the  configuration  of  the  network. 
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CAPTION 

27 

Figure  1.  The  real  and  imaginary  parts  of  in  dyne-cm/10 

Triangles  are  redrawn  from  (M;  Figure  6).  Squares  are  calculated 
by  matched  filtering  from  10  stations. 


The  Theoretical  Basis  for  the  Rapid  and  Accurate  Computation 
of  Normal  Mode  Eigenfrequencies  and  Eigenfunctions 


R.  P.  Buland  and  F.  Gilbert 


NOTE:  Computational  procedures  have  been  developed  for  the 

CDC  7600  at  the  Lawrence  Berkeley  Laboratory  of  the  University 
of  California. 


Introduction 


Recent  progress  in  the  application  of  normal  n.ode  theory  in 
seismology  has  made  it  necessary  to  develop  a practical  and  accurate 
scheme  for  the  computation  of  the  nigh  frequency  elastic-gravitational 
free  oscillations  of  a radially  symmetric,  non-rotating,  self  gravitating, 
perfectly  elastic  earth  model.  The  scheme  developed  is  based  on  the 
classical  variational  approach  and  on  recent  advances  in  the  solution 
of  the  algebraic  eigenvalue  problem  for  large  banded  systems  of  the 
form  (A  - AB)u  = £ . 

Raleigh-Ritz  Procedure 

By  Raleigh's  principle  the  eigenfunctions  of  an  earth  model 
are  extremal  solutions  of  the  energy  balance  equation : 

Q 3 O 3 

u f S'  (r)r“dr  - f (r )r “dr  = 0 (1) 

0 0 

2 

where  gj  is  the  kinetic  energy  density  per  unit  volume  as  a function 

of  radius,  */  the  corresponding  potential  energy  density,  u.  the 

angular  frequency,  and  a the  radius  of  the  earth  (Pekeris  and  Jarcsch, 
1958,  develop  eq.  1 in  terms  of  the  radial  scalars  of  an  eigenfunction). 
Therefore,  approximate  eigenfunctions  and  eigenvalues  may  be  computed 
from  eq.  1 by  a Raleigh-Ritz  procedure.  Let  us  represent  the  radial 
part  of  the  eigenfunction  s(r)  as  a linear  combination  of  N test 
functions  £^(r)  each  satisfying  the  boundary  conditions. 


2 


N 

s(r)  ~ s(r ) = 53  t.?.(r)  (2) 

i=l  1 1 

substituting  eq.  2 into  eq.  1 results  in  a matrix  eigenvalue  problem: 

A(u)b  = (u2T-V)b  = 0.  (3) 

with  the  following  properties: 

1)  Each  eigenvector  b^  represents  a projection  of  eigen- 
function s onto  the  space  spanned  by  the  C^'s  . 

2 

2)  u>  is  an  upper  bound  to  the  ’euared  eigenfrequency 
associated  with  eigenfunction  s . 

3)  Fach  sucessive  eigenvector  b_  represents  a higher  radial 
order  (overtone)  of  the  same  angular  order.  Each  angular 
order  is  represented  by  a different  matrix  equation. 

4)  If  the  error  between  s and  the  projection  of  s onto 

2 

the  C\'s  0(a)  then  the  error  between  u>  and  the 

2 

squared  eigenfrequency  assoc iatea  with  s is  0(a  ) . 

Test  Functions 

For  computational  speed  it  is  desirable  that  the  test  functions 
be  economical  to  compute,  that  the  number  of  test  functions,  N , be  as 
small  as  is  consistent  with  an  accurate  representation  of  s , and  that 
each  test  function  overlap  only  a small  number  of  other  test  functions, 
i.e.,the  matrices  T and  V be  banded  (Courant,  1943).  Further,  we 


demand  that  s be  continuous  and  desire  that  3 s(r)  be  continuous 

r 

also  except  at  first  order  discontinuities,  where  the  usual  boundary 
conditions  apply.  These  requirements  dictate  the  use  of  the  basis 
functions  for  piecewise  cubic  Hermite  interpolation  as  test  functions 
(Birkhcff  et  at.  , 1966). 


Computational  Philosophy 

In  order  not  to  discard  any  precision  unnecessarily  it  has 
been  found  desirable  to  adhere  rigidly  to  an  all  pervading  philosophy. 
As  is  common  practice  our  models  are  specified  at  discrete  points  of 
an  unequally  spaced  grid: 


v : 0 = Yl  < y2  < • • • < yM  = a . 

Each  of  the  model  parameters--  p(r)  , p(r)  , or  \(r)  --is  defined 
to  be  the  cubic  spline  interpolation  of  the  parameter  specified  on 
grid  v . Similarly,  the  eigenfunctions  are  calculated  on  grid: 

it  : 0 - x.  < x„<  • ' • < x,  = a . 

1 l • M 

Due  to  the  choice  of  test  functions  each  eigenfunction  must  be 
defined  as  the  Hermite  cubic  spline  interpolation  of  the  eigenfunction 
on  grid  ir  . 

It  has  been  found  necessary  to  perform  the  integrals  quite 
precisely.  Our  philosophy  has  teen  to  sacrifice  a little  speed  (since 
they  need  only  be  done  once  per  model)  and  do  the  integrals  exactly. 
Therefore,  all  integration  has  been  performed  interval  by  interval 
(xi_1  <xSxi,i=2j3,...j 


N)  by  means  of  a low  order 


^ 4 


Gauss-Legendre  quadrature.  The  integrals  are  thus  both  rapidly  computed 
and  numerically  stable  as  well  as  mathematically  exact.  It  has  also  been 
our  philosophy  to  use  no  less  than  6 nodes  (12  degrees  of  freedom  per 
radial  scalar)  per  radial  wavelength  at  the  highest  frequency  of  interest 
(G.  Frazier,  personal  communication). 


The  Detailed  Test  Functions 


Define  basis  functions 

px3  - 3x2  + 1 

0 < X < 1 

n0(x)  = n-2x3  ~ 3x2  + 1 

-1  < x < 0 

l 0 

otherwise 

jr  3 0 2 

fx  - 2x  + X 

0 S x < 1 

*Q(x)  = < X3  + 2x‘  + x 

-1  < x < 0 

L 0 

otherwise 

nQ  and  are  simply  designed  to  be  piecewise  cubic,  continuously 

differentiable  functions  with  properties  giver,  in  Table  1.  The  test 
functions  are  then  defined  by: 

fn  /»  -*i— \ x.  s X S x_ 

r°\xi+i  - v 1 1+1 

/ x - xi  ) x.  . s x < x. 

1 \xi  - xi-l / 1-1  x 


HjU)  = 


otherwise 


(5) 


t|/i(x)  = 


, S X s X.  , 
1+1 


( <xi,l  - ”i)  *0(^x7)  *i 

l*i  - -i-i)  M:  Vi  5 4 < xi 


0 


otherwise 


i 


Representation  of  the  Eigenfunctions 


Let  us  represent  a toroidal  free  oscillation  as: 


«£<:>*>  2 tnVr)  f x Vi  <e  •»>]  * 


. m. 

1 U)„t 
n l 


N 


nwp<r)  ~ [w^^n.ir)  + w^.2V(r)] 

i=l  * 1 11 


(2) 


(6) 


Eq.  6 yields  a matrix  system  with  2 degrees  of  freedom  per  ncce  or 
2N  x 2N  matrices  with  half  bandwidth  4 . Of  course,  for  a model 
with  zero  viscosity  in  the  outer  core, the  model  need  only  include  the 
mantle  for  this  case. 

Let  us  represent  spheroidal  free  oscillations  as: 


ni'E-.t)  = D?nVr>v>  ■ »>  * r.Yr)V°(9  , ,)]  e*n 


m. 

1 V 


„Ut(r)  - E [u/^np.)  YU.‘2\,<r)] 


(7) 


V.(r)  ~ E [v.(3)n.(r)  tv.l2)^(r)l 

1 x i l l J 


(2) 


coupled  with  the  perturbation  cf  the  gravitational  potentia] 


nY”(r.t)  = [fnPt(r)YX’(0,*)]  e 

N 


• nr 
x u.t 
n S, 


nPi(r)  ~ E [Pi(l)ni(r')  - p^2)  *.(r)] 


(2) 


(8) 


3.  .. 


jj 


3>6 


Eq,  7 alone  yields  a matrix  system  with  4 degrees  of  freedom  per  node  or 
4N  x 4N  matrices  with  half  bandwidth  8 . Eq.  7 coupled  with  eq. 6 yields 
a 6N  x 6N  matrix  system  with  half  bar.uwidth  12  . 

Boundary  Conditions 

Before  the  matrix  eigenvalue  problem  in  eq.  3 car.  yield  valid 
approximations  to  the  eigenfrequencies  and  eigenfunctions  of  an  earth 
model,  all  boundary  conditions  must  be  satisfied.  Rather  than  introduce 
constraint  equations  with  undetermined  Lagrange  multipliers  (which 
increases  the  size  of  the  matrices)  we  have  chosen  explicitly  to  match 
the  boundary  conditions  by  making  linear  combinations  of  the  rows  and 
columns  (and  reducing  the  rank)  of  the  matrices.  Therefore,  upon 
entrance  into  the  eigenvalue  procedure  all  boundary  conditions  are 
automatically  and  exactly  satisfied. 

Finding  the  Eirer.values  and  Eigenvectors 

By  using  the  Sturm  count  (number  of  eigenvalues  cf  eq.  3 

larger  than  u)  and  det  La(w)],  all  eigenvalues  in  a given  frequency 
band  may  be  quickly  found  to  macr.ir.e  precision  by  bisection  and  linear 
interpolation  (Martin  ana  Wilkinson,  1967;  Peters  and  Wilkinson,  1969). 
Even  though  each  calculation  of  the  determinant  and  Sturm  count  of 
eq.  3 requires  a new  decomposition  of  the  matrix  A(u.)  , it  is  possible 
to  take  advantage  of  the  bandwidth  of  the  system.  Without  the 
small  bandwidth  of  A(o-)  the  problem  would  be  hopelessly  expensive. 

Once  the  eigenvalue  is  known  the  eigenvector  may  be  found  by  inverse 
iteration  (Wilkinson,  1965,  p 628-629). 


J 


1 


A. 


If  *he  half  bandwidth  is  m , then  the  number  of  operations  needed 


to  find  each  eigenvalue  goes  as  Nm  . The  amount  of  computation  is 


independent  of  n,  l,  ui,  or  the  radial  turning  point  of  the  mode. 


Discussion 


The  unusual  feature  c:  this  method,  that  the  amount  of  computa- 


tion is  independent  of  the  mode,  : s due  to  our  ••xtrore  conservatism. 


No  advantage  can  be  taken  of  the  smoothness  of  lev,  frequency  modes  ci 


of  the  shallow  ray  equivalent  turning  point  of  the  high  angular  orcer. 


lov,  radia.  crdei  mantle  modes.  Cn  the  other  hand,  as  the  computation 


has  been  ~uie  practical,  the  luxury  of  such  conservatism  has  dene 


away  with  troublesome  special  cases.  Core  modes  and  Ctor.eley  modes 


are  as  easy  to  compute  as  any  other  type: 


Furthermore,  this  conserve t isr.  makes  possible  the  solution  of 


the  algebraic  eigenvalue  problem  with  the  following  advantages: 


1)  Highly  accurate  eigenfre juenc ies . 


2)  -orrect  computation  ci  nearly  coincident  * i/enr re  .uer.c  i ■ 


3)  All  modes  in  a frequency  _and  are  found. 


** ) prior  interred  tier,  about  eigenfrequencies  or  eigenfunctions 


is  required. 


The  eigenvalue  procedures  routinely  yield  eigenvalues  with 


a precision  cn  the  order  of  a part  in  10  . However,  the  accuracy 


3f  the  eigenvalue  ;3  0(h6)  , h - m*x(x.+1  - x.)  , and  the  accuracy 


of  the  eigenfunction  is  0(h)  (Eirkhcff,  et  at.,  1969).  For  our 


lbO  point  grid  ••  , we  evaluated  the  accuracy  emperically  by  comparing 


“D8 


eigenfrequencies  between  the  variational  program  and  a differential 
equation  program  (Alterman,  Jarosch  and  Pekeris,  1953).  The  differ- 
ential equation  program  was  written  with  10  decimal  place  constants 

0 

and  required  a precision  in  the  integration  of  cne  part  in  10 

The  eigenvalues  were  identical  to  a relative  acc_racy  of  no  less 

7 

than  one  part  in  10  . Since  these  methods  are  philosophically 

very  different,  and  algorithmically  completely  distinct,  this  ’esult 
indicates  to  us  that  both  programs  are  free  of  errors. 

As  a single  example  of  the  use  of  these  methods  we  present, 
in  Figure  1,  the  (w  , l)  diagram  of  spheroidal  modes  (a.  s 2m  -*5  sec 
1 S 150)  for  model  1Q56B  (Gilbert  and  Iziewonski,  la7E). 


! 


Caption 


Figure  1.  The  spheroidal  (w  , £)  diagram  for  Model  1066B 

(ui  ^ .14  sec.  , l < 150).  A continuous  line  joins 
the  points  for  each  fixed  value  of  the  radial  over- 


tone number . 


00- 0 
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Summary 

The  Poisson  sum  formula  is  used  to  transform  the  standing  wave 
representation  of  seismic  displacements  in  a radially  stratified  sphere 
into  a traveling  wave  representation.  The  terms  in  the  sum  can  be  given 
a physical  interpretation,  via  the  method  of  stationary  phase,  as  classi- 
cal wave  packets  making  successive  traversals  of  the  great  circle  through 
epicenter  and  receiver.  The  traveling  wave  representation  forms  the  basis 
for  the  retrieval  of  structural  parameters  and  source  mechanisms  from  the 
spectra  of  traveling  wave  groups. 
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The  representation  of  seismic  displacements  in  terms  of  traveling  waves 


1.  The  traveling  wave  representation 


The  representation  in  terms  of  normal  modes  used  by  Gilbert  and 
Dziewonski  (1975,  Sec.  2,  hereafter  referenced  as  M ) is  a standing 
wave  representation.  Ihe  sums  over  angular  order  can  easily  be  con- 
verted  into  integrals  cvar  wavenumber  from  which  a traveling  wave 
representation  can  be  obtained.  Such  a representation  is  desirable  for 
the  study  of  regional  structure,  for  the  separation  of  traveling  wave 
orbits  and  overtones,  and  for  the  study  of  the  moment  rate  tensor  of 
regional  events. 

Suppose  U(r,w)  is  one  component  of  the  displacement  spectrum 
for  a particular  radial  ' ^rtone  number,  n , (M,  [2.1.24,25]) 

CO 

U(r,u))  = £ a (r,«)  C#(w) 

a 1=0 

a^(r,u)  = A^(r)  • f(u >) 

and,  for  simplicity,  suppose  that  A (r)  depends  only  on  P (cos  0) 

* Jt 

aa(r,w)  = h(£+Js,  r,  u)  P^cos  0) 
h(i+is,  r,  u>)  = H^(r ) • f(w) 

A^r)  = H^(r)  Pt(cos  0) 


2) 


Then, 


£2 


There  are  two  steps  to  represent  U(r,u>)  in  terms  of  traveling 
waves.  The  Poisson  sum  formula  (Titchmarsh,  1943,  p.  60;  Nussenzweig, 
1965,  p.  27)  is  used  to  convert  the  sum  in  (3)  into  a series  of  integrals 

■ 00 

4)  u<r,w)  = / dA  £(-)k  h(  A ) p , (cos  0)  C(A,u>)  e'i2knA 

0 k=-«°  ^ 


and  zero  is  added  to  (4)  in  a suggestive  manner.  We  use  the  Legendre 
function  of  the  second  kind,  Q^_^(cos  6)  , and  rewrite  the  sum  in  (4) 

00  00 

5)  U(r,u)  = / dAh(A)  £ R (A, 6)  C(A,w) 

0 S=1  S 


where 
s - odd : 


RsU,e) 


6) 


s-  even: 


Rs(A,6) 


= (-)^s  1^2  [pA  l5(cos  0)  cos[(s-l)nA] 

♦ 7 Qa -Ig  ^cos  0)  sin[(s-l)irA]J 

s/2  r 

= I P^_^^cos  9)  cos(sttA) 

- I-  (cos  0)  sin(snA)J 


£3 


The  term  in  (6)  involving  Q^cos  0)  is  zero  for  s=l,  and  the  rest 
vanish  in  pairs  (2,3),  (4,5),  etc.,  so  that  the  terms  in  (5)  involving 
Q^cos  0)  sum  to  zero. 

The  reason  for  using  (5)  rather  than  (4)  is  that  we  can  interpret 
Rs(A,0)  in  (5)  and  (6)  as  representing  the  s--  passage,  or  orbit,  of  a 
traveling  wave  group.  This  interpretation  is  clarified  by  using  the 
pair  of  functions  (Nussenzweig,  1965,  p.  89;  Robin,  1958,  pp.  237-240) 

7)  Qa^2)(cos  b)  = % [Pa_15(cos  0)  ± i | Qa_^(cos  0)J 

and  rewriting  (6)  in  the  form 
s-  odd: 

»,<».•>  ■ (— {!'1W[q‘^(cos  e)e-i<s-1)rt.Q'2)(coS  6)e1(s-1)rt] 

8) 

s-  even: 

Rs(i,9)  = <->s/2  [q‘2’(cos  e)e-is*A  ♦ q‘^Ccos  B)eis”X] 

2.  A physical  interpretation 

To  understand  the  physical  interpretation  of  R (1,0)  , as 

s 

representing  the  s passage  of  a traveling  wave  group,  we  use  the 
method  of  stationary  phase  ana  the  asymptotic  approximations;,  valid 
for  Xt  » 1 and  e s 0 s tt  - e , 


1 


9) 


£ 4 


Qx-^2^Cos  0)  = (2,t*  sin  O)”*5  e^^0-11^4^  (1+0(A_1)  ) 

Also  we  consider  f(w)  to  be  a constant  in  (1)  and  (2).  The  actual 
behavior  of  f can  be  recovered  by  convolution  in  time  or  multiplica- 
tion in  frequency.  In  this  case  the  temporal  behavior  U(r,t)  corre- 
sponding to  (5)  is  (M, (2. l.ll)) 

10)  u(r,t)  = / dAh(A)  £ R U,e)  |cos(u( A )t ) e'°U)t  -l]  H(t) 

0 s=l  L J 

In  (10)  we  consider  the  transient  term  which  we  write  as 

ID  U(r,t)  = .>?„  / dAh(A)  f)  R (x,e)  eMA)t-a(A)t 

0 s=l  s 

Let  us  consider  s-  odd  in  (8)-(ll)  and  let  U(1’2)  denote 
the  contribution  to  U according  as  we  use  q(1^2^(cos  6)  in  (9). 
Then  the  integrand  for  IT1*2*  has  phase  $(1,2^ 

12)  ^1’2^  s w(A)t  + A [( s - 1 ) tt  + e]  ± tt/4 

According  to  the  method  of  stationary  phase  the  dominant  contributions 
come  from  those  values  of  A for  which  d$(1)/dA  = o or 

tdu(A)/dA  = (s-l)ir  + 0 = ^ i 0 

s 


13) 


If  we  define  a group  velocity  ys  = Ag/t  , in  units,  say,  of 
radians/s.,  then  (13)  requires  that 

14)  du(x)/dA  = A /t  = v > 0 

s Ts 

in  order  that  /1)  be  stationary.  For. toroidal  modes  the  variational 
formulation  (Meissner,  1926)  shows  that  du»(x)/dA  £ 0 , Thus  has 

at  least  one  stationary  value  for  tome  . for  spheroidal  modes  it 
is  almost  always  true  that  du(X)/dX>  0 . Some  contrary  examples  have 
been  given  by  Gilbert  (1967), 

The  interpretation  of  A in  (13)  is  that  it  is  the  total 
distance  traveled  by  the  wave  group.  The  s—  wave  group  has  traveled 
a distance  0 from  source  to  receiver  plus  an  additional  (s-l)/2 
great  Circles  in  time  t at  group  speed  ys  • This  is  exactly  the 
interpretation  given  to  the  classical  surface  waves  G.  , R , etc. 
Consequently,  s corresponds  exactly  to  the  orbital  index,  1,3,- 
for  such  a surface  wave.  In  addition,  the  s—  term  in  (11),  through 
the  factor  (-)  ^ in  (8),  is  out  of  phase  with  the  s ± 2n<*  term. 

This  i’ lustra tes  the  polar  phase  shift  (Brune,  Nafe  and  Alsop,  1961). 

The  s—  wave  group  has  passed  two  more  poles  than  the  s-2nd  wave 
group,  and  has  its  phase  advanced  tj/2  twice. 

Proceeding  in  a similar  manner  for  in  (12)  we  find  the 

stationarity  condition 

15)  du(X)/dX  = -A  /t  = -y  SO 

s s 


This  condition  is  never  met  for  toroidal  modes  and  only  rarely  for 

/ o \ 

spheroidal  modes.  Therefore,  U is  of  only  minor  importance,  and 
is  customarily  neglected. 

However,  in  (11)  we  could  have  used  exp(-iw(A)t)  in  place 
of  exp(iw( \)t)  . Then  U would  have  provided  the  dominant  con- 
tribution and  U(1)  would  have  been  negligable.  That  is,  the  dominant 
contribution  to  (11)  for  u>0  and  s-  odd  is  l/1*2*  according  as 
the  time  behavior  is  taken  to  he  exp(±ia»t)  . 

Considering  s-  odd  in  ( 8 ) -( 11 ) shows  that  Rs(A,0)  represents 
the  s traveling  wave  group.  The  group,  of  wavenumber  A and 

frequency  00(A)  , travels  a total  distance  A , given  by  (13),  in 

s 

time  t at  group  speed  y . The  important  term  in  R (A,0)  in  (8) 

s s 

has  superscript  (3,2)  according  as  the  time  behavior  is  exp(lio)t)  ; 
co  > 0 . 

Let  us  novf  consider  s-  even  in  (8)-(ll),  and  let  l/1’2*  be 
the  contribution  to  U(r,t)  according  as  we  use  Q^£2*(cos  0)  in  (9). 
The  integrand  for  l/1*2^  will  have  phase 


16)  *(1’2)  = u>( A)t  ± X [stt-0]  1 i»/4 

(2) 

The  phase  $ has  the  stationarity  condition 


17)  tdoj(x)/dA  = sir-0  = A >0 

s 

or 

18)  du(X)/dX  = a /t  = v i 0 

s ' s 


1 
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It  is  clear  that  the  s--  wave  group  travels  the  distance  2*-0 
from  source  to  receiver  plus  an  additional  (s-2)/2  great  circles  in 
time  t at  group  speed  Yg  . The  successive  values  of  s correspond 
exactly  to  the  orbital  index  2,  4,  • • • , f or  such  classical  surface 

waves  as  , etc.  The  polar  phase  shift  between  successive  orbits 

is  given  by  the  factor  (~)S^  in  (8). 

A consideration  of  leads  to  (15)  with  Ag  given  by  (17), 

and  U(1)  is  of  only  minor  importance.  The  rfiles  of  U(1)  and  U(2* 
are  reversed  if  the  integrand  in  (11)  is  replaced  by  its  complex  con- 
jugate. That  is,  the  important  term  in  Rs(X ,6 ) in  (8)  has  superscript 
(2,1)  according  as  the  time  behavior  is  exp(±iwt)  ; u>0 

The  foregoing  discussion,  based  on  the  method  of  stationary 
phase  and  the  work  of  Brune,  Nafe  and  Alsop  (1961),  illustrates  that 
(5)  is  the  exact  traveling  wave  representation  of  seismic  spectra. 

The  functions  Rs(X ,0 ) in  (8)  represent  traveling  wave  groups. 

3.  The  standard  form 

At  this  point  it  is  desirable  to  put  (5)  into  a form  similar 

to  (2.1,24)  and  (2.1.25)  of  M . Using  Legendre  functions  of  the  second 
kind 

19)  Qx^’2)(cos  «)  = h [‘"x.^cos  0)  ± i i Q^Ccos  0)J 

we  introduce  the  functions  Rm(A,0) 

s 9 


3-r  odd: 


.6 ) 

2C) 

s-  even: 


r"(X,0) 


()(s-l)/2 


(-) 


s/2 


[qa(,^(cos  e)e'l(s-1)"X  + Q'J^cos  0)  ei(s-1HAJ 
[q'^Ccos  9)  e-is"‘  t o’^cos  9)  eis’*] 


and  we  define 

21)  x"(At0) 


( -)m(  A/2TT 


(A24s)->5(A2-9/4)",s 


•r"(x,0) 


according  as  rn  = 0 , 1 ,2  . The  second  and  third  terns  in  braces  in 
(21)  can  be  well  approximated  by  A and  X*^  , respectively,  for  X>>1. 

We  now  define  As(A,r)  in  terms  of  (21)  exactly  as  A is  defined 
in  (2.1.30)  and  (2.1.31)  of  M in  terms  of  X*^(9)  . Note:  the  sign 

of  in  M(2.1.31)  should  be  negative. 

Thus,  we  can  write  the  traveling  wave  representation  generalizing 

(5),  as 


22)  U(r,w)  = / dX  £ a (X,r)  C(X,W);  a (X,r)  = A (A,r)  • f(u) 

0 s=l  S 9 S 

The  traveling  wave  representation  (22)  forms  the  basis  for  the 
retrieval  of  structural  parameters  and  source  mechanisms  from  the  fre- 
quency spectra  of  traveling  wave  groups.  It  also  forms  the  basis 


£ 9 


g 


for  an  investigation  of  body  waves  (Burridge,  1966;  Ansell,  1973), 
where  Ag(\,r)  is  decomposed  into  a series  of  terms,  each  of  which 
represents  a generalized  ray  (Brune,  1964;  Ben-Menahem,  1964). 

4.  Transition  to  plane  stratification 

The  representation  in  terms  of  traveling  waves  for  a plane 
stratified  medium  can  be  obtained  from  the  foregoing  expressions  by  a 
limiting  process.  Let  x = a0  where  a is  the  radius  of  the  sohere 
for  which  ve  have  the  traveling  wave  representation.  Let  x , the 
epicentral  distance,  be  fixed.  Let  ,\  = ka  , where  k is  the  (fixed) 
wave  number  on  the  surface.  Then  the  representation  for  a plane  strati- 
fied medium  is  found  as  the  limit  as  a-*"*°  of  the  representation  for  a 
sphere.  We  use  the  uniform  asymptotic  approximations  of  Szego  (1934) 
that  permit  us  to  replace  (9)  by  (23)  for  large  X 

23)  Q^2)(c°s  6)  = Ve/sin  e)H  hJ>2’1)(A0)  (itCU-1)  j ; 9<*/2 

in  terms  of  Hankel  functions. 

The  terms  for  s > 1 in  (5)  represent  wave  groups  that  arrive 

at  successively  later  times,  tg  , where  t is  the  earliest  group 

arrival  time  for  the  s—  group.  Since  lim  t = - , because  each 

a-**  s 

group  must  travel  a distance  of  at  least  tra  , we  shall  consider  only 
s = 1 in  (5) 

m 

24)  U(r ,<i>)  = / dXh(A)  R^A.0)  C(A,w) 


i 

i 

1 


1 


i 


1 


L.  v 


k. 
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Combining  (23)  and  (8)  gives 


Rj(M)  = *s(0/sin  e)*5  £h£2)(A0)  t 


= (e/sin  6 )”  Jq( X0 ) 


The  product  A0  in  (25)  is  the  same  as  kx 


26)  R^(A ,0 ) = (0/sin  0 P JQ(kx) 


As  0 -►  0 (a-*»)  (24)  becomes 


OD 

27)  U(r,o))  = a / dk  h(ka)  J (kx)  C(k,ui) 

f\  » 


A careful  analysis  of  h , including  the  normalizing  integral  (M,  (2.1.3)) 
shows  that  h is  0(a_1)  so  that  lim  of  (27)  ,is  finite. 


The  presence  of  the  cylindrical  Bessel  function  in  (27)  in  place 
of  the  spherical  harmonic,  the  Legendre  function,  in  (1)  and  (2), 
illustrates  the  transition  from  spherical  stratification  to  plane 


stratification. 
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